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ABSTRACT 

We introduce a collection of statistics appropriate for the study of spinorial quantities defined in three di- 
mensions, focussing on applications to cosmological weak gravitational lensing studies in 3D. In particular, 
we concentrate on power spectra associated with three- and four-point statistics, which have the advantage 
of compressin g a large number of t ypica lly very noisy modes into a convenient data set. It has been shown 



^ty 

previously by iMunshi & Heaverisl (|2009|) that, for non-Gaussianity studies in the microwave background. 



such compression can be lossless for certain purposes, so we expect the statistics we define here to capture 
the bulk of the cosmological information available in these higher-order statistics. We consider the effects 
of a sky mask and noise, and use Limber's approximation to show how, for high-frequency angular modes, 
confrontation of the statistics with theory can be achieved efficiently and accurately. We focus on scalar and 
spinorial fields including convergence, shear and flexion of 3D weak lensing, but many of the results apply 
for general spin fields. 



INTRODUCTION 
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this paper, we consider a set of new statistics designed to encapsulate much of the information content of third-order and higher statistics for spinorial fields 
de fined in three dimensions. In cosmological applications such higher-order statistics can be very noisy, and the dimensionality of the space may also lead to 
f^vfery large number of data to consider. Thus some form of data compression is attractive, but preferably in a way which does not reduce the cosmological 
,j^formation content inherent in the original statistics. In this paper, we build on the ideas originally presented in iMunshi & Heavens! dioO^ . where it was 
^own that one statistic, the power spectrum associated with the bispectrum, could be used very effectively to estimate non-gaussianity in the microwave 
. Jc^ckground radiation, as it was a lossless compression for this purpose, and also had the important added benefit of being able to provide evidence that a 
K»6n-gaussianity is primordial. In this paper, we extend the ideas to cover spin-weighted fields which are defined in three dimensions, with particular emphasis 
jqJi weak lensing fields convergence, shear and flexion. 

Ci ' Weak gravitational lensing of background source galaxies is caused by fluctuations in the intervening mass distribution. It manifests itself in a number 
of ways, most notably as distortions in their images. This effect arises due to the fluctuations of the gravitational potential and consequent deflection of light 

by gravity. 

Despite being a relatively young subject weak gravitational lensing i Munshi et alj|2008[) has made major progre s s within the last deca de, since the first 
measurements were published feacon. Refregier & Ellisll2000l ; IWittman et alll200(]| : iKaiser. Wilson & LuDpinoll2000l : IWaerbeke et "3120001) . There has been 
considerable progress in analytical modelling, technical specification and the control of systematics. By its dependence on the mass power spectrum at lower 
redshifts, weak lensing surveys play a complementary role to the studies based on large-scale galaxy surveys and Cosmic Microwave Background (CMB) 
observations. Ongoing and future weak lensing surveys such as the CFHT legacy survefl Pan-STARRsQ, the Dark Energy Survey, and further in the future, 
the Large Synoptic Survey Telescope^, WFIRST0and Euclid[f|will provide a wealth of information in terms of mapping the distribution of mass and energy 
in the universe. 

Owing to the lack of photometric redshift information the traditional approach to weak lensing has largely adopted a 2D approach, analysing correlations 
of the shapes of galaxy images on the sky only. However the availability of photometric redshifts a l lows a 3D weak lensing analysis, which was introduc ed by 



ot tne snapes ot galaxy images on tne sKy only. However tne availability or pnotometric redsnitts a l lows a iu weaK lensing analysis, wmcn was mtroauc ea by 
iHeavensI ^20031) . Later developments by various authors jHeavens et al.ll200ol : lHeavens et a]||2006l: Heavens. Kitching & Verdell2007 : Castro et al 20051) have 



shown that it can play a vital role in constraining the dark energy equation of state JHeavens et a |2006|) and the neutrino mass jKitching et al.l(2008h . This 
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has lead to recent progress in mo delling weak lensing observables in 3D extending results previously obtained in projection or using tomographic techniques 
l lMunshi, Coles & HeaveiislbOlOl) . 

Early results on analytical modelling typically assumed a small survey size and adopted a 2D approach that uses a flat-s ky formalism. This is relate d 
to the fact that first generation of surveys typically covered a small portion of the sky and lacked any redshift information l ljain. Seliak & White! I2OO0I) . 
Indeed such analytical modelling was very successful in p redicting lower-order statistical properties of weak lensing convergence and shear very accurately 
l lMunshi & Jainl200ll ; lMunshill2CI0(]|: lMunshi & JainllOOO). These results depends on analytical modelling of underlying de nsity pert urbations using perturba- 
tive a nd empirical methods. jMunshi & Jain 200ll : IValageasll200dl : lMunshi & Valageail2005l : IValageas. Barber. & Munshi 2004; Va lageas. Munshi & Barbed 
I2OO5I) . A tomo graphic step was next advocated to tighten the cosmological constraints. The tomo graphic studi es typically divides the sources into a few 
redshift slices i IHuII 19991 ; iTakada & Whit3l2003l : lTbkada & Jaia,2004. ; .Massey et al 2007 ; Schrabback et alj|2009l) . These slices are then ana lyzed essentially 
using a two-dimensional approach but including the correlation between different redshift slice s. A notab l e exce ption to the 2D analysis was IStebbinsI ( [l99a) 
who developed an all-sky formalism for weak lensing surveys. The techniques developed in IStebbinsI ( 1 19961) re ly on a tensoria l for malism, whe r eas w e 
will be using an equivalent treat ment based on spin weight spherical harmonics. Extending previous studies by iHeavenj ( l2003h and ICastro et a 3 l l2005t) . 
iMunshi. Heavens & Coles! bOlOl) extended the all-sky formalism to 3D to take into account the photometric redshift information, as well as extending to 
higher-order statistics. However they focused on the convergence field which, being a spin-0 field, is relatively easier to analyze. The main motivation behind 
this work is to extend previous results to arbitrary spinorial fields such as shear and their derivatives flexion. 

Weak lensing at small angular scales probes the nonlinear regime of gravitational clustering, and the extra modes there can lift degeneracies about 
backgroun d cosmology present in s tudies involving the po wer spectrum alone see, e.g.,[Bernardeau, Van Waerbeke & Mellier ( 1997); Jain & Seljak ( 1993); 
IHuH ( I1999I) ; !Schneider et al ( !2002h ; !Takada & JaiiT ( !2003!) . The nonlinear regime is characterized by gravity-induced non-Gaussianity, and detailed studies 
that employ the Fisher matrix formalism have already demonstrated the potential of using higher-order non-Gaussianity information to lift cosmological de- 
generacies. Higher-order studies are also important in evaluating the vari ance of lower-order statistics, e.g. a proper knowledge of the trispectrum is essential 
for computing the error bars in the power spectrum jTakada & J ain' 2009). The modelling of higher-order statistics typically involves either perturbative tech- 
niques or ernpirical modelling of the underlying m at ter clusterin g (Fry 1984; Schae f fer 1984 ; Bernardeau & Sch aeffer 1992; Szapud i & S zalay 1993, 1 9971 ; 
!Munshi et all 1 9991 ; !Munshi, Coles & Melotll999iJlbl ; !Munshi, Melott & ColeJl999l ; !Coles. Melott & Munshil!l999l ; !Munshi & Coles!!200d!2002ll2003l) . Us- 
ing such pre scriptions and their exterisions , studies involving non-Gaus sianity, have also been performed in projection (2D) as well as using tomographic 
information ( !Hull999! ; lmada & JairJ!2004!2003l ; [Semboloni et a]|2008h with remarkable success. 

Studies involving h igher-order correlation functions have been performed using observational data jBernardeau. Van Waerbeke & Mellieil !l997! ; 
!Bemardeau, Mellier & Van Waerbek e 2002; Pen et a l 20031 ; [jarvis, Bernstein & Jain!!2004!) . Most of these studies involve one-point moments (cumulants) 
which collapse the entire correlation function into a single number. Mode-by-mode estimates of higher-ord er correlation functions or multispectra though far 
more interesting is difficult given the low signal-to-noise of current observational data. Current studies bv !Munshi & Heavens! l [2009l) defined power spectra 
associate d with each multispectrum that u ses an intermediate option in data compression. While initially this concept was applied to CMB studies, recent 
work by jMunshi. Coles & Heavens!! 20 id) extended this concept to weak l ensing. This initial work focused o n convergence k. Being a spin-0 (scalar) object, 
the analysis of convergence statistics is relatively simple. In their analysis [Munshi. Heavens & Coles! j2010l) used the similar statistics for shear and flexion 
fields but in projection (2D). Th main motivation for the present study is to use the full 3D information (available from photometric redshift surveys) in 
analyzing the non-Gaussianity not only in the convergence field but also in shear and flexion. This is particularly interesting as current photometric redshift 
surveys with good image quality will provide a wealth of data for the analysis of weak lensing which can be used to probe cosmological information. For 
our study, we combine well-motivated ansatze in modelling the gravitational clustering with the Limber approximation. The results that we derive here are 
generic and will be useful in other areas of cosmology where integration along line of sight is involved. To keep the results simpler we will ignore the fact, 
that in a realistic survey, the average density of sources will decline with distance, and the distance estimated from photometry will also include error, but 
these are evidently important ingredients in a practical implementation of these statistics. 

The expressions for higher-order multispectra generically include multidimensional integrals involving multiple spherical Bessel functions. We will be 
using Limber approximation to simplify these results. We will show that, at each order, we can reduce the dimensionality of these integrals to unity by using 
Limber approximation. This will simplify the numerical evaluations of such integrals considerably. 

Thi s paper is arranged as follows. In §2 we discuss the basic formalism of 3D weak lensing. The formalism presented here is a generaliza- 
tion of i Munshi, Coles & Heavens! !201Q!) and jMunshi. Heavens"& Coles! !2010!) and can analyze higher-or der statistics of spinorial fields in 3D. In 
jMunshi, C oles & Heavens 2010) results were derived for higher-order statistics for the convergence and in ( iMunshi. Heavens & Colesl!2010l) the focus 
was on higher-order statistics of spinorial objects but in projection (2D). The notations for 3D harmonic decomposition, which will be used in the following 
sections are also introduced here. In §3 we introduce the models describing higher-order clustering of underlying matter which are then used to construct 
models for the bispectrum and trispectrum in the nonlinear regime. The results obtained are generic and can describe higher-order statistics of weak lensing 
convergence, shear and flexions. In §4 we focus on power spectra associated with higher-order multispectra. Results presented in this section correspond 
to both all-sky and patch-sky coverage. In §5 we focus on error analysis and derive results for scatter (or variance) of various estimators in the presence of 
observational noise and mask. Finally §6 is devoted to discussion of the results. Though we have mainly focused on weak-lensing, the general formalism 
developed in the paper will have wider applicability. We will use the Hierarchical ansatz to model clustering of underlying mass distribution, but the treatment 
can also be adopted in the context of more elaborate scenarios of clustering e.g. halo model. 



2 NOTATION 



This s ection is devoted to introduc ing the basic no tation and formalism of 3D weak lensin g . We w ill follow the n otation used in Munshi. Coles & Heavensl 
l!20IOl) which is based mainly o n [Heavens! l !2003l) and further developed by !Castro et all l l2005h . The results of ICastro et all ll2005l) were generalize d by 
Munshi. Coles & Heavens 1 20 IC ) to take into account higher-order correlations. The aim of this paper is to extend both Munshi. Coles & Heavens! 1 I2OIOI) and 



Munshi. Heavens & Coles! ( 120101) to the analysis of shear using full 3D information. 
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2.1 A Tale of Two Potentials 

Linking tlie 3D lensing potential and the 3D gravitational potential $ is crucial in connecting lensing observables to theory. In this subsection we will 
consider the harmonic decomposition of 3D scalar (spin-0) fields which is a step towards making this connection because examples of such fields include the 
scalar potentials and the convergence field k that we encounter in weak lensing. The harmonic decomposition is most naturally done using eigenfunctions that 
can be constructed using ordinary spherical harmonics and spherical Bessel functions. In the next subsection we will generalize them to the case of spinorial 
fields. 

The statistics of shear and convergence can be expressed in a natural way through their relation to ^{r,9, tp) the 3D gravitational potential at a 3D 
position r, 9, (p, and i^(r) the lensing potential. The density contrast (5(r) is directly related to the potential through the Poisson equation. This allows us to 
link directly the statistics of the weak lensing observables to the underlying statistics of the mass distribution, and hence to cosmological parameters. The 
radial distance r{t) is related to the Hubble expansion parameter H{t) = a/a by r{z) — c dz' / H(z'). The Hubble parameter is sensitive to the contents 
of the Univer se thereby making weak lensing a useful probe to study dark energy. The line of sight integral relating the two potentials can be written as 
teaisej|l992h : 

0(r)^0(r,n) = 4 r^^'^K(r,r')$(r',n); F^{r,r') ^ J^^^^—^ (1) 
Jo [./K(r)/K(r')l 

The Born approximation was used to derive the above expression iBemardeau. Van Waerbeke & Melliejl 19971 ; [Schneider et"a il l2002| j Waerbeke et a il l2002l) 
The lensing potential 0(r) = 0(r, Cl) has a radial dependence and is a 3D quantity. In our notation r = r{t) is the comoving distance to the source 
whose observed light was emitted at a given instance of time t. The observer is situated at the origin. The function FK(r, r') depends on the background 
cosmology, through the function /k(''); fK{r) = sinr, r, si nhr for a closed (K = 1), flat (K — 0) or open (K = —1) universes respectively. Our 
convention for the Fourier transform for the 3D fields is as in iMunshi. Coles & HeavensI bOlOl) . The eigenfunctions of the Laplacian operator in flat space 
when expressed in spherical coordinates turn out to be a product of spherical Bessel functions ji{kr) in the radial direction and the spherical harmonics on 

the surface of a unit sphere i.e. Yim(fl) = Yim,{9, </>). The eigenfunctions Zkim{r) = k ji{kr) Yim{^) are associated with eigenvalues —k^. In general 
the radial eigenfunctions are ultra-spherical Bessel functions, but they can be approximated by spherical Bessel functions when the curvature is small. The 
eigendecomposition and its inverse transformation can be expressed as: 



= / d'v^{r)ZM^{r)- (2) 



and 



oo m = l p 



(3) 



The specific choice of eigenfunctions allows allows us easily to expr ess the coefficients of expansion of the convergence (or shear) in terms of the 
expansion of the density field through the Poisson equation ( lHeavensll2003h A"I>(r) — 3QmHQS{r) /2a. In the harmonic domain this can be expressed as 
$imik;r) — ASim{k;r)/a{r) k with A = —3QmHo/2. Here, <I>im(fc) is the spherical harmonic decomposition of "I'(r), and similarly for ^(r). In our 
notations, a{z) — 1/(1 + 2) is the scale factor at redshift z, Qm is the total matter density at 2: = 0, and Ho is the Hubble constant today. Simik; r) is the 
eigendecomposition of 5{v). When appearing after the semi-colon, the r dependence (e.g. of (fc; r)) is really an expression of the time-dependence of the 
potentials, which translates to a dependence on r, as r depends on look-back time. Using these deco mpositions, the har monic decomposition of the lensing 
potential cj>im{k) and the 3D gravitational potential $im(fc, r) are related by the following expression jCastro et afcOOSh : 



Ak roc nr 

4>im.{k) ^ — 5- / dk'k' / r'^drji{kr) / dr' FK{r,r')ji{k' r')(^im{k' ■,r') 

TTC^ Jo Jo Jo 



(4) 



The basis functions for the harmonic decomposition of the spinorial fields such as flexion and shear will involve spin-weight spherical harmonics which we 
will introduce next. The 3D power spectra for the gravitational potentials $ and the lensing potential (f> are defined through the following expressions: 

(<l>i„(fc)<l>,,„,(fc')) =Cf*(fc)5iD(fc + fc')4'5^m'; (<^;™(fc)0,'™'(fc')> =Cf^fc,fc')5/r<Sw. (5) 



2.2 3D Eigendecomposition of Spinorial Functions 

In this subsection we will introduce the generic spin-weight functions and their eigendecomposition. Specific cases that are of interest here include shear and 
flexions. This will generalize the spin-0 results discussed above for the convergence field. We can expand the fields such as shear 7±(r), flexions J-{jc), Q{r) 
in 3D basis functions that are constructed out of spin-weight spherical harmonics sYim(Q) on the celestial sphere and spherical Bessel functions ji(kr) in 
the radial direction. Expansion in such bases provides a very simple relationship between harmonic coefficients of the shear, flexion and convergence on 
the one hand, and the lensing potentials on the other. Moreover, spherical coordinates are the natural choice for eigendecomposition as this provides a clear 
separation in terms of radial modes and the modes on the surface of the sky, and the ubiquitous presence of a sky mask induces mixing of modes only on 
the surface of the sphere, and the use of photometric redshift estimates only introduces error in the radial direction without altering the angular position. 
The choice of eigenfunction is also motivated by the Poisson equation which relates the 3D potential $(r) to the density distribution 5{r), whose statistical 
property we will model to predict the statistics of shear, convergence or flexion. Extending the definition of spin-0 eigenfunctions Zkim{r) we will denote 

the spin-s eigenfunctions as aZkim.{r); which is defined as: sZkimjr) = ^f^kj i{kr) sYimj^) ■ The spin-weight spherical harmonics are defined in terms 
of D-matrices jVarshalovich. Moskalev &Khersonskiilll988l :|p enrose &RindleJ[l984. 19861) . They satisfy a orthogonality relationship similar to ordinary 
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Figure 1. Left panel compares the power spectrum C; for the convergence k for two different source redshift Zs = 1 and Zs = 0.5. Right panel shows the power spectrum 
C; for the convergence k, electric part of the shear E, flexions and 5 as a function of I. The flexion power spectra Cf ,Cf are normalized by for display. The cosmology 
assumed is ACDM and all sources are assumed to be at the same redshift of Zs = l.The A CDM background cosmology that we have considered is described by the following 
set of parameters: Cm = 0.3, Ua = 0.7, T = 0.21, h = 0.7 and as = 0.90. 



spherical harmonics. Spin-weight harmonics with the same or different spin indices are orthogonal on the surface of sky. This generalizes the orthogonality 

relationship of ordinary spin-zero spherical harmonics; sYim{^) ~ \J'^^D^-s,m{^^ 0) 0), J dtl aYim{^) siYv m' {^)d(l — SuiSmm' Sss' ■ 

Alternative expansion schemes are indeed possible such as using tensor spherical harmonics but they are perhaps more difficult to wo rk with. It is wort h 
noting here that the formalism of spin-harmonics is extensively used in studies involving Cosmic Microwave Background Polarization l lBunn et af2003l) . 
The forward and inverse transform of an arbitrary spin function s/(r) from real space to harmonic space links it with its harmonic components s.ftm{k) that 
can be expressed as: s/(r) = dkJ2\Zo' $I^m=-it/im(fc)]sZfc;„i(r) and afim{k) = / d'V[s/(r)]s^fci„(r). The orthogonality relationship satisfied 
by the 3D spherical basis functions [sZli^{r)] depends on the orthogonality of spin- weight harmonics sYim{^) and that of the spherical Bessel functions 
ji{kr). It generalizes a similar relation for the scalar harmonics. For arbitrary spinorial fileds with spins s and s' it reads: J v[a Z'li^{v)] [^1 Zyii^, (r)] = 
Soik — k')5iiiS^^i Sggi . The inverse transforms are used to define the harmonic components of generic spinorial fields ?7(r) and ??*(r). The results that 
we will derive in our later sections are expressed most naturally in the harmonic domain using these components 2Vim{k) and -2ilim{k) which can be 
expressed as: 2'nim{k) = / d^r r]{r)2Zkim.{r) and -2'nim{k) = / d^r r?*(r)2^feim(r). It is indeed possible to work with 2'nim{k) and -2riim{k) as well 
as the harmonics Eim{k) and Bim.{k) that can be constructed from them. Though they contain the same information, the Electric or E and Magnetic B 
modes provides a rotationally invariant description in full sky. The expansion coefficients Eim. has a parity ( — 1)' while Bim has a parity ( — 1)'+^. The clear 
separation of modes with different parity gives a clear mathematical advantage in the case of weak lensing, as It can be shown that, at first order, weak lensing 
from gravitational clustering can only generate E modes, whereas systematics are mostly responsible for the generation of any B mode contribution. 

The explicit expressions for the electric Eim{k) and magnetic Bim.{k) components, constructed from these harmonic transforms are: Eim{k) — 
— |[2r7im(fc) + -2T]im{k)] ; Bim{k) = |[2?7im(fc) - -2?7im(fc)]; and ±2riim{k) = —[Elmik) ± iBim{k)]. The individual components of the field 77(r), 
r]\ (r) and 7)2 (r) are expressed in term s of eigenfunctions Z+.kim { v) tha t can be constructed from linear combinations of Z±2,kim {v) introduced before. The 
formalism used here is very similar to lMunshi. Heavens & ColesI ( l2010h . The emphasis here however is not just on 2D decomposition on the surface of the 
celestial sphere but rather on a 3D decomposition which relies on the photometric redshift to estimate radial distance. 

It is worth mentioning here that unique decomposition of a function into modes E and B mode on the celestial sphere is possible only with complete sky 
coverage. In the presence of a boundary, which is often the case owing to the presence of masks, the decomposition is ambiguous. For the case of weak lensing 
shear these equations can be specialized further by ignoring the magnetic contribution which is zero for shear generated pur ely by gravitational lensing in the 
absence of any systematics. Indeed, higher-order lensing corrections can generate lensing B mode too l lCoorav & H j|2002h but are sub-dominant. 



2.3 Harmonic decomposition of Convergence, Shear and Flexion 

The results derived in previous section can directly be applied to the case of flexions, shear and convergence. Most of the generic results are applicable to the 
analysis of shear 7 if we specialize the field with a spin— 2 object and identify with 3D shear 7. Complex shear 7 constructed from its individual components 
7±(r) — 7 i (r) zb ^72 (r) acts as a spin-2 object and can be expressed in terms of the lensing potential 4> using spin-derivatives (see[Munshi, Heavens & ColesI 
1 120101) and lCastroet"Sl bOOSi) for more discussion on spin-derivatives) which are used to construct spinorial fields with different spin-weights. The lensing 
potential plays the role of the generic scalar field introduced earlier to express arbitrary spin functions. We will use the same symbol <j) for both. We will use 
the generalized symbol sF for general spin fields which will include products of shear fields as well as higher derivative spin objects such as flexions. In our 
cuffent notation 2r = 7 and _2r = 7*: 2r(r) = 7(r) = i99[(/!>] and _2r(r) = 7*(r) = i96[0*]. In general the scalar potential ^(r) will have both 
electric ^^(r) and magnetic (t>B{r) components: (f){r) = 4>e{v) + i(l}B{r). 

The individual shear components 71 (r) and 72 (r) and the convergence K(r) can be expressed in terms of a complex lensing potential (j){r) — (f>E{jc) + 
icj)B{r). As pointed out before the magnetic part of the potential (fiBijc) will take contribution mainly from systematics and the electric part corresponds 
largely to pure lensing contribution 71 (r) = i( 9 9 + 9 9 )4>{^) \ 72 (r) = —\{^ 3 — 9 3 )0(r) and k{y) = ^(9 9 + 9 9 )(/>(r). 

Derivatives of the shears are higher-spin objects. Using these derivatives quantities such flexions are constructed, which are also often used in the 
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Figure 2. The diagonal components of the reduced convergence bispectrum bm are plotted as function of I. The bispectrum is assumed to have a hierarchical form. The 
hierarchical amplitude Q3 is set to unity. The source redshift is set to unity Zs = 1. The plots for the flexions and Q are normalized by l'^ for display. 

contex t of weak lensing studies jGoldberg & NataraianI ( l2002h : lGoldberg & Bacoi] ( l2005l) : lBacon et al.l ( l2006h : lBacon & Goldberg ( l2005h : ISchneider cfe Eil 
The two flexions that are most commonly used are also known as the first flexion T (spin — 1) and Q which is also known as the second flexion (spin 
3). These two flexions in combination can specify distortion beyond what is described by shear. The flexions can be used to describe weak "arciness" in image s 
of lensed galaxies and their relationship with the shapelet formalism is well documented l lRefregieil2003l ; lBernstein & Jarvisl2002l ; lRefre gier & Bacon 20031). 
The flexions J-{r) and 0(r) are both used in the literature mainly for individual halo profiles and also for the study of substructures l lBacon et al. 20061) . 
We are mainly interested however in higher-order statistics of these objects for generic underlying cosmological clustering. This is done by linking the 
3D harmonic decompositions of the flexions to that of the lensing potential 4>{r): J-{r) = i(609 + 09S + 999) (j>{r) and C/(r) = i 9 9 9 4>{r). 
Flexions have been used primarily to measure the galaxy-galaxy lensing to probe the galaxy halo density profiles. Their cosmological use will depend on an 
accurate understanding of gravitational clustering at small angular scales. 

In Fourier space the harmonics of 7(r) and 7*(r) can be expressed in terms of the harmonic coefficients of $B(r) and $s(r) denoted by Ei„^{k) 
and Bi,n{k) respectively: ±2Tim{k) = —[Eim{k) ± iBim.(k)]. Analogously, the harmonics of T and Q denoted by J-im{k) and Ghn{k) can also be 
expressed in terms of the (f>i,n{k). In the absence of B-modes the harmonics of the shear components are directly related to the harmonic component of the 
Electric field Eim- The harmonic transforms of the shear components and convergence can also be expressed in terms of the lensing potential as follows: 

^im{k) = -'-^^ci>l^{k)■ Ei^{k) = _iy/l|±||^,„(fc); Ti,,,{k) = ll^'^{l + lf/^ [il^ + il-2)(l>i^{k); = iy|±|[0,„(fc). 

These harmonic expressions can be used to reconstitute the real space spinorial fields: ±2r(r) — J kdk X];^o \J ('-2)! 4'im{k) ±2Zkim{r, fi); 

J"(r) = / kdk J^t^Q Ylm=-i ^im{k) -iZkini{r) and C/(r) = J kdk J2uLo Qim{k) aZkim{r). These results derived above are useful in finking 

the statistics of weak lensing fields K(r) 7(r), J-{r) and G(r) with those of the underlying density field <5(r) responsible for generation of the lensing potential 
^(r) with the help of Eq.®. 

The theoretical modelling of the underlying mass distribution that we employ in our study is based on the hierarchical ansatz. The hierarchical ansatz is 
more suited to model gravitational clustering a smaller scales, which makes it particularly suitable for modelling the flexion statistics which put more weight 
on smaller scales. A comment about noise contribution due to intrinsic flexions of source galaxies is in order. While it is relatively easy to model the intrinsic 
ellipticity of source galaxies, detailed modelling of intrinsic flexion of source galaxies is much more complicated and depends heavily on modelling of galaxy 
shapes beyond the simplest description. This uncertainty is also expected to increase with survey depth. 

In addition to shear, convergence and flexion which are used in weak lensing studies, we can also consider a generic scalar tracer field >!' in our study. 
Such fields can represent a suitable large scale tracers which are sometimes used for cross-correlation studies or studies involving weak lensing magnification. 

The statistics of shear and flexions can be best related to that of convergence with certain / dependent multiplicative factors that we will call form factors. 
In later sections Ff will denote the form factor associated with a generic spin-weight field F. So the form factor for the shear 7± will be denoted by F^^ . 



3 3D WEAK LENSING STATISTICS : POWER SPECTRUM AND BEYOND 

For the study of non-Gaussianity we need to go beyond the study of power spectra. In this section we will present results for 2-, 3- and 4-point statistics, 
and show the relation between observables and theory. The various multispectra involve multidimensional integrals, which we simplify by employing various 
levels of approximations involving the high I behaviour of j;(a;). 

3.1 Power spectrum 

We will start by deriving the power spectrum (fci , k2) for the 3D weak lensing fields. Our derivation of the 3D convergence power spectrum is based on 
expressing its harmonic coefficients Ki,n{k) in terms of the 3D density field S with the help of Poisson's equation and using the definition of the convergence 
field in terms of the projected lensing potential <j) gives: 
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2kA, 



K,m{k) = —1(1 + 1) dk'k' r^drjiikr) dr' Fy^[r,r') ji{k' r') 



i.5im.(k';r') 



k'^a{r') 



(6) 



We will use a shorthand notation Ii{ki, k) (defined below) useful for si mplification of our results. We will approximate the cross-spectra at two different 
epoch using the approximation P**(A:,r,r') = JP'^'^lk. r]P'^'^'(k. r') jCastro et allboOSi) . Use of this approximation leads to separation of respective 
integrals. As we will see below the use of the extended Limber approximation, which is valid at high /, implies that dominant contribution will come from 
single time slices r = r' and this approxim ation is not detrimental to any of the final results which are quite generic. The 3D power spectrum can be expressed 
in terms of Ii(ki, k) as lCastro et all ( 1200^ : 

oo n r 

Ii{ki,k) = h drr^jiihr) dr' Fi,{r,r') ji{kr')./P^^ky) 
Jo Jo 

16 /'°° 1 ' ' 

Cf'^{ki,k2) = k^I,{ki,k)Ii{k2,k)dk; Cr(fci,fc2) = -/'(/ + l)'cf'^(fci,fc2); Cr {ki,k2)^FfFf Cr{ki,k2). (7) 
Jo 4 

Clearly the above expression is quite generic and contains all the weak lensing information at the second-order level. This expression is however quite 
cumbersome for any numerical implementation as it involves three-dimensional integral which are quite demanding computationally. We will be using 
extended Limber approximation valid at high I to simplify the above expression. Using this approximation we can reduce the integrals to one-dimensional 
integrals. In any case we will quote the generic result that is valid without any approximation. Notice that the following approximation is also independent of 
the factorization of the power spectrum introduced before. 



cr(fci,fc2) = 



16 



k^dkli{ki, k)Ii{k2, k) 



16 



kik2 drarijiikiT-a) dr'aFK{ra,r'a) drbrlji{kirb) dr'i,FK{rb,r'i,) k'^dk P'^'^{k,r'a,r'i,) ji{kr'a)ji{kr'i,) 



(8) 



We will next use the Limber approximation Ea.l lA2t to simplify the k integral which produces a Smir'a — J"f,) function. Integrating out r'^ with the help of 
the delta function and renaming the dummy variable r'a to r' we can finally write: 

Ci'^{ki,k2) ^ -ki k2 I rldriji(kiri) j rldr2ji(k2r2)T>f'' (ri,r2); 

12 









'Di'"'{ri,r2) = ^ [ r'^-^^ FKiri,r')FK{r2,r')Ps ( ^■,r' ) ; r„i„ = mm(ri , ra). (9) 
c Jo v ) V 



Use of the Limber approximation projects multi-time coiTelators to a single time correlator. Going one step further. If we use the high I approximation to the 
spherical Bessel function Eq. dASt to reduce the dimensionality of the above integrals involving the spherical Bessel functions ji, we arrive at the following 
simpler approximate equation. Use of Eq.l |A3t allows us to replace ri and r2 in terms of fci, fea and /. 



2,2 



cr{ki,k2) = [Ac-'] 



21 + 1 



21 + 1 



2ki 



21 + 1 



2k2 



,2 dr' 

^ 21 



21 + 1 
2ki 



-,r 



Fk 



21 + 1 



2k2 



(10) 



We can define a statistic E(fci, fca) which will include all available information from individual harmonics, as a function of fci, fca: 



fca) = Y.^21 + l)Cr (fci, fca); E'"'(ri, ra) = ^(2? + l)Cr (n, ra). (11) 
I I 

We have ignored angular smoothing in our derivation. Typically observations will involve a smoothing filter. Tophat and compensated filters are the ones 
that are most commonly used that can be incorporated in . As pointed out before the above equation is derived using very general arguments. It is valid at 
high I as the derivation is based only on high I approximation to the spherical Bessel function ji [x). Nevertheless the derivation of a 3D skew spectrum has 
wider applicability in cosmology. The technique can be applied to compute 3D power spectrum in other context (e.g. integrated Sachs-Wolfe effect or Kinetic 

Sunyaev-Zeldovich effect). Detailed analysis for such cases wi ll be presented elsewhere. 

We have used Limber approximation to simplify results ILo Verde & AfshordU ( 1200 Sh . It was pointed out by ILo Verde &Afshordil l l2008h that using I 
instead of / + i, as is often done in the literature, spoils the accuracy of Limber approximation to 0{j). In general the error in Limber approximation will 
scale as ©(ji). A series expansion of spherical Bessel function can also be performed to construct the next to leading order terms which further improve the 
accuracy of Limber approximation. 

Weak lensing not only induces correlations among ellipticities of background galaxies (shear), it also introduces local modification in the number 
density of source galaxies (also know as weak lensing magnification). The weak lensing magnification ^ is directly linked to the convergence k. While we 
have focused on shear 7± in this paper we plan to extend the results to magnification in a related publication. 

While the results derived above are valid for all-sky surveys, observations invariably will introduce mask. We will next analyze the case of 3D power 
spectrum estimation in the presence of a general mask. The results that we derive will have general applicability and will be valid for near all-sky coverage. 
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Figure 3. The skew spectrum c'^'"''^ defined in Eq. )3.U for the convergence ft is plotted as a function of I for two different source redshifts Zs = 1 (upper curve) and Zs = 0.5 
(lower curve). A hierarchical form for the bispectrum was assumed and the hierarchical amplitude is set to unity. As ACDM background cosmology is assumed. We have not 
incorporated any smoothing window - calculations are done by introducing a sharp cutoff at Imax = 2000. 



3.1.1 The Effect of an angular sky Mask 

We start by 3D decomposition of an arbitrary spinorial field sF in the presence of a sky mask w{Cl), whicli is equal to or 1 in simple cases. The decomposition 
into radial harmonics, using spherical Bessel functions jiikr), can be performed independent of the mask. The harmonic decomposition on the surface of 
the celestial sphere involves spin-weight spherical harmonics sYim,{d). The following expression relates the masked observed harmonics [sTwlimik) and 
he unmasked s^imi^)- We will leave the spinorial field gV arbitrary and derive the expression of the cross-correlation power spectrum, in the presence of a 
mask, with another arbitrary spinorial field s/r'(n). The results are generic and do not depend on any specific assumption that is used to model the all-sky 
power spectrum itself. 



E 



kr'dkji{kr) / dn[sT[r)w{Q.)][sYi'^[Q.)] 



dk k ji{kr)sTi>m'{r) 



§ E E i-^y^"" Jdr Jdk'k k' ji{kr)ji,{k'r)sTi,^,{k')wi^^Jvi^i 
+„(2r + l)5 ^ 



I' m' li^ mfj 

E E (-1) 



(2Z + 1)2 



21 + 1 



21' + 1 



lUlamaJl'lal 



I' la I 

s -s 



I' la I 

S -S 

I' la 

m' rria 



I' la 
m' IJla 



I 



(12) 
(13) 



The masked cross-spectrum Cf^ (k, k') involving T(Cl) and r'(f2) can be described in terms of the all-sky cross-spectra Cf" (k, k') and a mode-mixing 
matrix G that describes the effect of mode-mode coupling resulting from the presence of the mask. The mode-mixing matrix depends on the spin-weight of 
the respective fields s and s' and also depends on the power spectrum of the mask wi . 



C7\k,k') = -l—Y^[:T]Uk)[s~r%..{k') = -l—Y^[gTwU{k)l.T'wU{k') « ^QrC^' 



2/ + 1 ^ 21 + 1 



21' + 1 ' 21' + 1 



1 ^p (2r + l)^ (I la I' 



I la I' 

s' -s' 



(14) 
(15) 



The radial direction remains unaffected by the mask which introduces mode mixing only on the surface of the celestial sphere hence the matrix M is 
independent of radial wave number k. For a given pair of radial wavenumbers k, k' the all-sky cross-spectra Cf^ {k, k') can in principle be recovered 
by inverting the above expression Ea.lll5t. In general there will be contribution from noise which can be from intrinsic ellipticity or flexion distribution of 
galaxies in case of shear or flexion. Such contributions need to be subtracted to make any estimation unbiased, and there may be standard issues with inversion 
which may require regularisation. 

To recover the power spectrum of E and B modes of a spin ±2 fields, commonly used in the context of analysis of CMB polarization analysis, e.g. 
as in lBrown. Castro & Tavlod jlOOSh , we have to express angular harmonics of sF =+2 F and ^/F' —-2 F in terms of their Electric (E) and magnetic (B) 
components 2F;„i = Ei^ + iBim. Then using Eq.OSj we can relate the cut-sky power spectra Ci = Ylm, EimE^^ and C; = jITT BimBim. 
in terms of their all-sky counterparts Cf^ and Cf^. However Eq.l|15|l generalizes such results to generic spin functions with arbitrary spin-weights. For 
generic spin weight functions the cut-sky and all-sky relations are: 
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CEE f^EE ^EE I y^iEB^BB ^BB ^BE^EE . ^—iBB^BB 

A sum over repeated indices are assumed in eacli of tliese equations. Tiie matrices GfiF , GfiF and GfiF are defined througli tiie following expressions: 



It is interesting to notice here that instead of Ci{ki, ^2) if we study E(fci, k2) — X]i(2' + l)Ci(fci, k2) they will have exactly similar mixing properties 
as the ordinary 2D fields, modulo the remapping of the radial harmonics, as the usual 3D power spectrum C; when a mask is applied, i.e. E'^'^(fci, ^2) = 

Mill E ( 21,^1 ^1 ' 2/'+! ^2); where Ga' = '''^21+1) • ''"'^'^ property will be valid not just as the level of power spectrum but also for skew- and kurt 
spectra as well as for multispectra of arbitrary order. 

We have plotted the projected power spectrum for convergence Cf" in Fig. (D as a function of / (left panel). Two different redshifts were considered 
Zs = 0.5 and Zs = 1.0. The ACDM background cosmology that we will be using throughout this paper is characterised by the following set of parameters: 
fim = 0.3, Qa = 0.7, F = 0.21, h = 0.7 and erg — 0.90. The power spectra associated with other spinorial fields (shear and flexions) are also plotted (right 
panel). The 3D power spectrum C'^''''^ is plotted in Fig.|4]for the same background cosmology. We plot C'^'^''^{k, k) for three different choice of k values as 
function of / (left panel) as well as "^'"(fc, k) for three selection of / values as a function of the radial wave number k. 




3.2 Bispectrum 

The power spectrum carries the bulk of the information in any cosmological observations. However often a set of degenerate cosmological scenarios can 
lead to a very similar power spectrum. Analysing higher-order correlation functions can lift this degeneracy to some extent. The non-Gaussianity used can be 
either due to primordial or secondary effects, and in the case of weak lensing the main source of non-Gaussianity comes from gravitational instability. Note 
that a non-zero bispectrum signifies the lowest-order departure from gaussianity, and its detection is generally easier than higher-order multispectra. 

To make contact with the observables we use the fact that the convergence can be related directly to the 3D density field. We will start by linking the 3D 
convergence bispectrum B and the 3D density bispectrum expressed in harmonic coordinates. In the next section we will express the bispectrum in spherical 
coordinate in terms of the bispectrum in rectangular coordinates and use some well-motivated approximations to simplify the results. 

Statistical isotropy requires that 

(Kiimi(fci;ri)Ki2„2(fc2;r2)Ki3m3(fc3;r3)) = ^ B^ihhC'i'ri) (18) 

and using Eq.l|6j we can write 

B^^i.(k,r.) ^ A^C.C2C. (^) r f r <^r.rl,iMr'M^^^^^ ^^-(-.-) >< 

/ TT I dr2rlji2{k2r'2)ji^(k2r2) [ 4^ -Fk (r-2 , ) / ^ / drlrzji^{k'ir'^)ji^{k3r-i) I -^^fj-FK{r-A,r'3)Bf^i^i^{k'i\r'i)\ 

Jo '^2 Jo JO '^('''2) Jo '^3 ^0 JO "■Vs) 

A = U{U + 1) ~ it (19) 

The bispectrum Bi^i'^i^{ki\ri) can now expressed in terms of the underlying matter bispectrum . The above relation mixes modes only in the radial 
directions r, and on the surface of the sky there is no mixing of angular harmonics if there is no sky mask. While expressing the density harmonics in terms 
of the 3D potential harmonics, we pick up additional scale factor a{ri) and wavenumber ki dependence in the denominator. 

We have so far ignored the presence of noise. Indeed because of the limited number of galaxies available it may not be possible to probe individual 
modes of the bispectrum at high signal-to-noise ratio. In later sections we will be able to address issues related to optimum combinations of individual modes 
which may be better suited for observational studies. 

The convergence bispectrum can be written in terms of the density bispectrum as follows, using the Limber approximation to simplify the results: 



poo poo poc ^ 

I3'!'ll':^i^{ki;ri) ^ H1H2H3 rldriji^{kiri) rldr2ji2{k2r2) I rldr3ji^{k3rz)ll^]^i^{ri,r2,r-j,); Hi = A 
Jo JO Jo 

^l%i,{r,,r2,r,)^Si,i,i.Mi.H-Si,i,i, T""" r^drB* (^1, li, ii; r, r, r) 7^l(r)7^2(r)7^3(r■); n.{r) = ^^^^ 

Jo \r r r J a(r) 

{2h + I){2l2 + l){2h + 1) f h h h 



Shi.i.-]/ ^ J. (20) 

To derive this result we have used the extended Limber approximation Eq.( IA2b to simplify the k'^ integrals. The integral here extends to the overlapping 
region i.e. rmin = min{r\,r2, r-i). In particular we can use the gravity-induced bispectrum here, or include others, such as a primordial bispectrum. 
It is possible to simplify further the above expression using Eq. dAlt : 
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(3) 



2h + 1 2/2 + 1 2/3 + 1 



2fci 



2fc2 



2fe 



(21) 



We will use a generic hierarchical ansatz to model the matter correlation hierarchy. Such ansatze constructs multi-point correlation functions from the products 
of lower-order correlation functions. In the Fourier domain this will lead to the construction of multispectra from products of ordinary power spectra. Such 
models have been tested against simulations and are routinely used both for projected galaxy surveys and for 2D weak lensing surveys. At the level of 
bispectrum we have: B{ki, k2, k-i, r) — Qs\P(ki,r)P(k 2,r) + P( k2, r)P(k-i, r) + P{ki, r)P{k3,r)]. More detailed mo delling will make Q3 a f unction 
of the wave vector triplet (fci , ^2 , fc3) e.g. in the halo model lCoorav & Seth ( :200Z) or in Hyper Extended Perturbation Theory jScoccimarro et alll998h . While 
the expression derived above is for the convergence field k, it can be used to construct the other bispectra involving shear or flexion: 



Results involving flexion can be constructed replacing the form factor F,^* s with the ones for the flexions i.e. or Pf defined accordingly. The radial 
dependence of convergence, shear or flexion harmonics are the same. 

Mode coupling is introduced by the presence of an observational mask. It is not possible to deconvolve the effect of a mask while analyzing the 
bispectrum from a realistic survey, as the inversion is typically unstable. However later we will introduce a power spectrum associated with the bispectrum 
(the skew spectrum), which can be computed from realistic data in the presence of mask and deconvolution can be done in a way very similar to the estimation 
of power spectrum discussed previously. 

If we assume the intrinsic shear and flexion of source galaxies to be distributed according to a Gaussian distribution, then they do not contribute to the 
estimated bispectrum, but this is, of course, an assumption. However even in this case the scatter or variance in estimation does get a contribution from such 
a source of noise. 

It is important to realize the results derived above are generic. They do not depend on the specific model used as an example (hierarchical ansatz). If we 
replace the underlying bispectrum with a primordial bispectrum of a specific type (e.g. local) we can still use the formalism developed here to compute various 
relevant statistics, we will introduce later, e.g. the skew spectrum. Later we will also introduce an optimized estimator for the skew spectrum. This estimator 
is not only optimized to detect any specific type of non-Gaussianity but it can also give an estimate of leakage from a specific source of non-Gaussianity (e.g. 
gravity-induced) while estimating another (primordial). 



3.3 Trispectrum 

The trispectrum or, alternatively, the four-point correlation function can provide an important sanity check jKamionkowski. Smith & Heavensll201Clh to 
validate lower-order detection of non-Gaussianity based solely on the bispectrum. The trispectrum, being a four-point correlation function is generally harder 
to probe compared to the bispectrum. However, in cases where the bispectrum vanishes due to symmetry considerations, the trispectrum is the lowest probe to 
study gravity-induced non-Gaussianity. In addition, while the study of trispectrum is interesting in itself it is also important in the proper characterization of 
the error in the power spectrum. As in the case of bispectrum we will start by modelling the bispectrum of the convergence field which can then be generalised 
to model trispectra associated with various spinorial fields. These results will eventually be used to model two different power spectra associated with the 
trispectrum (the hurt spectra). 

The convergence trispectrum T^^i^ {L, ki \ ri) is the four-point correlation function in the harmonic domain and can be expressed as 

rr4V4(L,fc,;rO. (23) 

LM ^ ' ^ ' 

The vectors l\,l2, h, h represents the sides of a quadrilateral and L is the length of the diagonal. The matrices as before are the Wigner 3j symbols. The 
symbols are only non-zero when they satisfy several conditions; which are |/i — /2I < L < /i + I2, j/a — /4I < 1/ < /a + h; h + I2 + L= even, I3 + U + L 
= even and mi + m2 = M as well as ma + 7714 = —M. In our notation for the trispectrum, TI^^^ {ki,ri\ L), the indices (ki-ri) encode their dependence 
on various Fourier modes of the density harmonics in the radial direction, used in their construction. No summation will assumed over these variables unless 
explicitly specified. To model the trispectrum we will relate it to the underlying trispectrum of the density distribution T^^l^ (fc^; r^). 

-Thh^i, ^ Air r r r f f ^''A H '^^'^ H ^ 2- (i.' '^ H p ( M 

T,X'{h;n)^ACrC2C3C,[-^^y--[-^^^ — | d™,,(/c,r,)| ^^-(n.rO 

r-dk^ /°°dr47|j,,(fclrl) r Jf±-Fj,ir,,ri) ^T^V^fc^rO. (24) 

Jo '^4 JO JO "l'"4j 

The above expression is a direct consequence of Eq.(|4ll. To make any further progress we need to consider a specific form for the matter trispecrum. There 
are two generic prescriptions for treating a gravity-induced trispectrum. The halo model is one, and has been developed over the last several years and is very 
popular for modelling the correlation hierarchy of the underlying mass-distribution. Another is that perturbative descriptions can also provide a reasonable 
description of the onset of non-linearity at comparatively larger length scales. The hierarchical ansatz on the other hand describes the matter correlation 
hierarchy in the highly non-linear regime on smaller length scales. It builds up the higher-order correlation hierarchy from the two-point correlation function. 
All possible diagrams that connect points at which the correlation function is being constructed are considered. These diagrams are attributed various 
amplitudes according to their topology. Different diagrams with same topologies are associated with same amplitude. At the level of the four-point correlation 
function there are only two different topologies star and snake. We will denote the corresponding amplitudes by Rt and Ra and consider each of these 
contributions separately next. 
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3. 3. 1 Star Diagrams 

The star diagrams are easier to handle because topologically they consist of a single vertex and lack any internal momentum that needs to be integrated out. In 
generic hierarchical scenarios a star diagram appears at each order in the hierarchy. In the case of projected 2D analysis, it has bee n found that replacing all dia- 
grams with the same numb er of star diagrams often is sufficient to reproduce all one-point statistical features of convergence maps l lMunshi. Valageas & Barbeii 
l2004l : lBarber. Munshi & V alageas 2004; Valageas, Munshi & Barber 2005). However the results presented here are generic and includes both contributions. 

The derivation of the star contribution to the convergence trispectrum follows a similar technique as the bispectrum. The first step is to express the density 
convergence in terms of the triplets of matter power spectra . Next, using the expression Eq.Q we can relate the star contribution of convergence trispectra to 
that of underlying mass distribution: 

^T}^}^ ^L,^-p-rt^ = kik2k3k4, J dririji^{kiri) . . . J dr4r4j^ (fc4r4) j/^'';*^;^^ (ri, r2, rs, r4)star, where 

Ji%i^i^{ki;r,)stiir = RbSi^i^LSi^i^L / r^drji^(kir) ... ji^{k4r) {P{ki;ri)P{k2;r2)P{k3- rs) + cyc.perm.} . (25) 



The above equation is for the stellar contribution to the trispectrum of the underlying density distribution. The additional three terms can be recovered by 
cyclic permutation of the ki variables. We can use this result next to express the convergence trispectrum. The simplification relies on the use of the extended 
Limber approximati on to simplify the kj integrals. H ierarchical ansatz and Limber approximations are both known to be valid at small scales which justifies 
their combined use l lMunshi. Coles & Heavensll2010l) . 



'^'^hu I ^ H1H2H3H4, f rjdriji^{kiri) ■ ■ ■ f r|dr 4^, Jfc4r4)X;**'^,^;Jr-i, r2, rs, r4)Btar where 

V /star ■'0 Jo 

1-\fi^i^i^{ri,r2,r3,rA)stB,, = RbSi^i^hSi^uL r'^drTZi{r) . . .7^4(r) |p {~'^^ ^ iv'^^ ^ iv'^^ ^ cyc.perm.j . (26) 

In generic hierarchical scenarios the trispectrum is a cubic combination of underlying matter power spectra, just as the convergence bispectrum is a quadratic 
function of matter power spectra. We will focus on a specific model, the hierarchical ansatz, to model the underlying density distribution. However as it was 
pointed out that similar techniques can also be applied in the context of more elaborate halo model prescription. 



3.3.2 Snake Diagrams 

The analysis of the snake diagrams is more difficult than that of the star diagrams. This is related to the fact that while higher-order star diagrams at each 
order are straightforward generalizations of the lo wer-order star diagrams, the simke diagrams are however constructed using two different lower-order star 
diagrams. The following expression was derived in lMunshi. Coles & Heaveri3 ( l2010h which relates the star contribution to convergence trispectra. In addition 
to various form factors, the following expression depends on the cu bic product of the underlying m atter power spectrum. The relative importance of star 
and snake differs in different models of the hierarchical ansatz, e.g. ISzapudi & Sz"aiavl l ll99l [l997l) attributes equal weighting to both diagrams, whereas 
iBemardeau & Schaeffeil l [l992l) constructs higher-order amplitudes from the lower-order star amplitudes. The amplitude of the new star diagra m at each order 
is left arbitrary which can be fixed using simulations or through the use of hyper-extended-perturbation theory jScoccimarro &Friemanlll999l) which predicts 
one-point cumulants at each order. We note in passing that ordinary perturbation theory predicts a hierarchy similar to hierarchical ansatz at the onset of 
gravitational clustering. However the amplitudes in this regime for various topologies are different ( iFrviil984) . 



/ / ■ \ f°° 

= H1H2H3H4 ridnji^ikin) ■ ■ ■ rldr4.ju{k4.r4)ll^}^i^i^{ri,r2,r3,r4.)sna.kc 

V /snake "'0 Jo 

^l^l^i^i^iri, r2, r3, r4)sm,kc = RaSi^i^LSi^i^L J r'drTZi(r) . . .TZ4{r) i^P (^^,r^ P (^^^,7-^ P (J^,r^ +cyc.perm.|. (27) 

The eye. perm, here represents a total of three other terms. These terms can be obtained by rearranging ks {h — >■ I2), {I3 — > I4) and (/i I2, 13 I4). 
The other terms that can obtained by considering two additional pairings by considering the exchanges {I2 I3) and {I2 -t- I4). These will lead us to Q'^'iJ 
and Qjgj*. The total number of snake terms considering three distinct pairings and permutations within each pairings is twelve. The snake contribution to 
trispectrum can be written as: 



^ / snake ^ / snake ^/ v J \ / snake 



L" 



The matrices in curly braces are the 6j symbols i lEdmonds|[T968l) . The differences in snake and star contributions are also apparent in various choices of 
permutations of (/i, hjhjh; L) associated with individual terms. 

The derivations outlined both for bispectrum and trispectrum are quite generic and depend only on the use of extended Limber approximation, and the 
analysis can be generalized to higher-order multispectra. 
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3D Power Spectrum 




Figure 4. The left panel shows 3D power spectram 1(1 + l)fc^C;'"' (fc, k) as a function of k for three different values of /. Where as the right panel shows 1(1 + l)A:^Cj'^'^ (fc, k) 
as a function of k for fixed value of I. A LCDM background cosmology was assumed with Qm = 0.3, Ha = 0.7, h = 0.7, erg = 0.8, Us = 1.0, w = —I; we use an Eisenstein 
& Hu (1996) linear power spectrum and the Smith et al. (2003) non-linear correction. The sharp cutoff in these plots reflects survey depth through the Bessel inequality. Here k 
is displayed in hMpc~^. 



We will use these results to construct the power spectra associated with the trispectrum, or kurt spectra. Construction of trispectra for generic spin-weight 
functions can be achieved by using the form factors in a manner similar to the one adopted for the bispectrum Ea.i22t. 




The two different kurt spectra that we will construct next will provide independent probes of the underlying mass trispectra, and can play a valuable role in 
checking any cross-contamination from systematics. 



4 POWER SPECTRA ASSOCIATED WITH MULTISPECTRA 

The multispectrum may encode a great deal of information, but there is certain amount of degeneracy involved in it. Owing to the low signal-to -noise 
associated with the estimation of multispectra, it is impossible to estimate multispectra mode by mode. Estimation of multispectra is also hampered by 
their complicated response to the survey mask and complicated noise characteristics. Recent studies have shown that degenerate sets of power spectra can 
be constructed from multispectra at a given order. These compress some of the available information in the multispectra and can be computed from the 
observational data with relativ ely ease. We will construct the p ower spectra associated with bispectrum and trispectrum in this section. These power spectra 
were studied in some detail in ( iMunshi. Heavens & Colesl201(]|) in projection (2D). We generalize these results here to 3D. First we will obtain generic results, 
applicable irrespective of detailed analytical modelling of underlying multispectra.Next we will further specialize these results for the models outlined above, 
and ee use extended Limber approximation to simplify the generic results. 

4.1 skew spectrum 

4.1.1 All-Sky Results 

We will start by expanding the product of two generic spinorial fields with associated spin indices s,s' [srgir']{r) in 3D in terms of their individual 
harmonics. The product of two spinorial fields with spins s and s' is a spinorial field of spin s + s'. The harmonic expansion therefore will have to be in terms 
of spin-weight spherical harmonics with spin index s + s'. In addition to expanding on the surface of the celestial sphere, we will also consider the expansion 
in the radial direction using spherical Bessel functions. 
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Figure 5. The right panel shows 3D power spectrum l{l + l)k'^C^'^''' {k, k) as a function of fc for three different values of Z. Where as the left panel shows /(i + l)fc|Cj'^'^''° (k,k) 
as a function of I for fixed value of k. Three different k values were chosen as depicted. A ACDMbackground cosmology was assumed with Qm = 0.3, Ha = 0.7, h = 0.7, 
(T8 = 0.8, Us = 1.0; we use an Eisenstein & Hu (1996) linear power spectrum and the Smith et al. (2003) non-linear correction. A hierarchical ansatz was assumed for the 
bispectrum with the amplitude fixed at Q = Q3 = 1. See text for more details. Here the wave vector k is displayed in units of hMpc~ ^ . 



We have used the Gaunt integral to express integrals involving spherical harmonics in terms of the Wigner 3j symbols at the last step. The matrices describe 
the Wigner 3j symbols. To construct the skew spectrum, next we contract the multipole [sT{r)s'r' {r')]tm with the multipole associated with a third spinorial 
3D field s"^'i'm{k3)- The resulting skew spectrum or power spectrum associated with the bispectrum depends on the three radial wave numbers (fci, k2, fca), 
in addition to the spinorial indices associated with the respective 3D fields: 



m 

The bispectrum is defined by a triangular configuration in multipole space with lengths of side {li,l2,l). The power spectrum constructed above, will 
essentially capture information, through a summation, with all triangular configurations with one of these sides kept fixed at length /. If we now use the 
following expression for the 3D bispectrum Bf^iji^ (fci, k2, ks) introduced earlier, 



Sfi^a {ki,k2,k3;ri) = ^ {sri^rni{ki;ri),,r'i^^^{k2;r2)s"Ti^m3ik3;rs))c ( 
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we can express the two-to-one power spectrum as 
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Although the above relation is general and contains all the information regarding the evaluation of the skew spectrum for arbitrary spin-weight functions, the 
presence of multidimensional integrals makes it difficult to implement computationally. 

We will simplify the results by using the extended Limber approximation to reduce the integrals involving fci and k2. First we express the B in terms of 
the underlying matter bispectrum B using the expression Ea.J20b. The extended Limber approximation, Eq. l lA2t . collapses the ki integrals to one-dimensional 
Dirac delta functions. These delta functions are next used to reduce the rj and r2 integrals. These simplifications which are known to be valid at high I. Finally 
we are left with two integrals involving r and 7-3 which can be further simplified by Eq. l lA3t . The final expression contains only a one-dimensional integral 
and directly relates the skew spectrum C]^^ (fc, fc') to the matter bispectrum: 
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Finally if we use the expression for the bispectrum derived using the same approximation, we can write: 
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The bispectrum B for convergence is defined in Eq.l|19ll. The generic symbol B^^iJ^ denotes the bispectrum for arbitrary spinorial fields as defined in 



Three-dimensional Statistics of Weak Lensing Shear and Flexion 1 3 



Ea.ll32t. Individual cases can be obtained by expressing B^^i^f in terms of the convergence bispectrum using the form factors introduced in Eq.l l22t . A 
numerical implementation of this result may not be easy due to associated computational cost. 

The diagonal elements of the reduced bispectrum bm are plotted as a function of I in Fig.|2] We have plotted the projected skew spectrum for convergence 
Qi..K,h Fig.[3l Two different redshifts were considered Zs — 0.5 and Zs — 1.0. The ACDM background cosmology that we have considered are described 
by the following set of parameters: Q,m — 0.3, Qa = 0.7, F — 0.21, h = 0.7 and erg — 0.90. We have taken a hierarchical form for the matter bispectrum 
and the amplitude Q3 is set to unity. Changing the amplitude Q3 only changes the overall normalization of the curve. Individual values of Cf^ at a given 
I depend on the modelling of the bispectrum for the entire range of / values being considered. Accurate modelling of the window function is therefore will be 
an important ingredient in such calculation. For our study we have considered a sharp cutoff in the multipole space at /max ~ 2000. The 3D skew spectrum 
Qh.ix,K plotted in Fig.|5]for the same background cosmology. We plot C^'"''^{k, k) for three different choice of k values as function of / (left panel) as well 
as Ci'^''^ {k, k) for three selection of I values as a function of the radial wave number k. 



4. 1.2 The Effect of a Sky Mask 

A mask on the sky is defined through a generic function w{Q) — X];™ wimYi„T,{i^) on the surface of the sky. The mask can be a simple zero and one step 
function signifying masked and observed part of the sky or it can also be a more complex apodizing function with specific weights attach to the different 
parts of the sky. We will compute the skew spectrum in the presence of mask and express it in terms of the skew spectrum in the absence of any mask. This 
will allow us to eventually design an estimator which can estimate the unbiased skew spectrum from the real data in the presence of a mask and noise. We 
will consider Gaussian noise, which means that the estimation of bispectrum is not affected as the noise has vanishing skew spectrum; however the scatter 
associated with the estimator will change in the presence of noise. We will denote the masked power spectrum by Ci. We will see that convolved or masked 
skew spectrum is a linear sum of individual all-sky C;s. This related to the fact that the use of mask introduces mode-mode correlation. A matrix inversion 
with suitable binning can produce the all-sky (7;s from the masked Cts. 

In our derivation we start by expanding the masked product field [sr(r)3/r'(r)TO(f2)] in 3D and contract the resulting harmonics with the masked 
harmonics of another 3D spinorial field [suT" (Cl)'w(Cl)]l^. The 3D harmonic expansion of the masked product field can be carried out using harmonics of 
spin index s + s' . Consequently the presence of a scalar spin-0 mask does not alter the spin of the product field which remains the same as the unmasked 
product field. 

Repeated application of Gaunt's integral allows us to write the product harmonics in terms of the individual harmonics: 
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The Fourier-Bessel decomposition of the quadratic field along the radial direction can be used to relate the 3D harmonics of the product field with that of 
individual constituent fields. This will eventually allow us to express the 3D skew spectrum in terms of the 3D bispectrum. The 3D skew spectrum presented 
here is a generalizati on of our previous work and incorporates full 3D information. This can be viewed also as a natural generalization of the 3D power 
spectrum presented in lCastro et all ( l2005l) . Higher-order counterparts at the level of trispectrum will be presented in the next section. 




.r(r),,r'(r)™(n)],™(fc) = 

3/2 



(:^) ^ Jhht'Jt'iaii-'^)'^'' J drkr^ ji(kr) ^ wi^^^ J dki ji{rki) ki sTi^^^ir) J dfe j; (rfe) s'Tj^^j (r) 
s s' —{s + s') J y mi m2 — m' J y (s + s') — (s + s') J y m' ma —m 



(37) 

Construction of the masked skew spectrum follows the same approach. As we pointed out before the 3D skew spectrum has a radial wave number k 
dependence built in it. If we integrate the radial dependance out we can recover the usual projected skew spectrum. Expressing the skew spectrum in terms of 
the multipoles of squared field and individual fields (suitably masked), 



1 ' 

{k, k') = [.r(f7),.r'(n)™(t7)],,„(fc)[,.r"(J7)™(f7)]L(fc, k'). (38) 

m—~l 

For the 3D harmonic decomposition [giiT" {il,)w{Cl)]*^{k) of individual fields we use the expression derived previously in Ea.lll3t. Using these expressions 
for the multipoles in the presence of a mask from Ea.J37t and Ea.lll3t and in Eq.(l38t we can express the skew spectrum constructed from the masked 
multipoles or the pseudo skew spectrum in terms of a coupling matrix and the all-sky skew spectrum in 3D and a mode-coupling matrix. The mode-mode 
coupling matrix involves only the mask power spectra wi = X^m "''"i^;?n' depends on the spinorial indices of the respective 3D fields: 

Guy'-r^pj^(^^^ + ^){sis' -isis')){l" s,s',s" eO,l,±2,S; (39) 

The spin indices s,s',s" take values for k, ±2 for 7±, —1 for F and 3 for G. After tedious but straightforward algebra along the line described in 
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iMunshi. Heavens & ColesI ( l201(]t) we can show that the pseudo-CfS expressed as a linear combination of all-sky power spectra can now be expressed using 
the following relationship: 



p,YT',T",, ■r-^^s',a"^rr',r" / 2/ + 1 21 + 1 ,\ 

{k,k)^}_^G„ C, ^^^k,^^k). (40) 

The deconvolution process to recover the all-sky C^,^ (ri, r2) from its masked counterpart involves inversion of the coupling matrix G which depends not 
only I indices but also on spinorial indices associated with the fields. In case of small-sky coverage which might be the case for present generation of surveys, 
for inversion of the matrix will typically involve binning to avoid any possible numerically-singular matrices. 

Here it is worth pointing out that the various spinorial fields that we can use to construct individual skew spectrum do probe the same underlying matter 
bispectrum. This can be used as a helpful diagnositic to probe possible spurious effects of mask and noise. 

The power spe ctrum {k, k') reported here is an extension of similar power spectrum introduced in lMunshi & HeavensI l l2009l) for CMB studies. 

Lat er it was used in Munshi. Coles & Hea vens (201Q;) to probe the convergence skew spectrum and was also generalized for projected spin-skew spectrum 
in (iMunshi. Heavens & Colesll2010h . In~this study we present skew spectrum for spinorial fields using a complete 3D analysis. The power spectrum is 
specified by multipole indices on the surface of the celestial sphere as well as radial harmonics along the line of sight direction. Results are generic for 
fields with arbitrary spin and can be used to probe shear, flexion or convergence maps. Similar statistics in the coordinate space has been reported before. 
ISemardeau. van Waerbeke & Mellieil ( l2003l) studied {'y^{Q,)^{Q,')) which directly deals with shear maps as opposed to convergence maps. This statistics 
along with a similar but simpler version which uses {h^{&,)k.{Q,')) was also studied. The perturbation theory was employed to model the underlying mass 
distribution and it used a flat sky approximation to simplify their calculations. A complementary statistic ([7(^1) ■ 7(02)17) was also considered which relies 
on more detailed modelling of the bispectrum. These statistics were used by ISernardeau. Mellier & Van Waerbeka (l2002h later to detect non-Gaussianity 
from the VIRMOS-DESCART Lensing Survey. 

Our results presented here deal with power spectrum associated with the higher-order multispectra, are derived using generic all-sky treatment, and can 
also handle decomposition into Electric and Magnetic components in a much more straightforward manner. The results presented here are not only applicable 
to shear or convergence but are also applicable for higher-order spinorials such as Flexions. To increase the signal-to-noise of the estimates it is customary to 
often sum all possible mode of C^^ {k, k) in to a single number which is called skewness. 



Sr''""' {k, k') = ^(2Z + l)Cr''^" {k, k') (41) 
I 

For a concrete expression for ^3 (k, k') we need to replace C]^^ {k, k') by the expression derived in Ea.ll35t or Eq.ll33t. 
To make connection with the real space statistics we can use the following relation: 
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and a similar relation can be derived for the skewness ^3^^ (r, r') and S^^ (fc, k') defined above. One point statistics such as S3{k, k') only contain 
radial information as all spherical harmonics are already summed over. 

If we make the approximation of replacing the spherical Bessel function with a delta-function form we can write 
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(43) 



An equivalent expression is valid for S^^ {k, k'). In both Eq.ll42b and Ea.ll43t we use the same notations to define the harmonic space C 
and real space C^^ {r, r') power spectra. However their functional dependence on their arguments are different. It is interesting to notice that in Fourier 
domain the C^^ (fc, k') probes specific length scales. 



4.2 The Kurt spectrum 

The four-point correlation function, or its harmonic counterpart the trispectrum, has been extensively studied in the literature . This contains the information 
about the non-Gaussianity beyond the lowest level i [HuI|I999I : lOkamoto & Hull2002l) . 

For the case of weak lensing studies clearly the gravity-i nduced non-Gaussianity is the m ain motivation. Studies in trispectrum analysis h ave also been 
pursue d using various other probes e.g. using 21 -cm surveys l lCoorav. Li & Melchiorrill2008h or more extensively in several CMB studies; see lBartolo et ^ 
( l2004h for a review. However these studies typically probe the trispectrum induced by primordial non-Gaussianity. It is important to note that at the level 
of four-point, the Gaussian part of the signal as well as the noise both carry a non-zero (unconnected) trispectrum. This degrades the signal-to-noise for 
various estimators and clearly needs to be subtracted out before an unbiased comparison with the theoretical predictions can be made. It is obvious that the 
detection of the trispectrum from noisy data is far more nontrivial than the estimation of the bispectrum. We provide analytical expressions here mainly for 
completeness and to show that the generic results can be obtained based on very simple hierarchical ansatz. 

Previous studies have mainly concentrated on one-point estimators which collapse the data to a single number - known as the kurtosis. We extend 
studies involving kurtosis {[sV ^'T' s'"T"']{(l,r)) to its two-point counterparts: {[,T ^,V']{Cl,r)[snV" ^,„V"']{Cl' ,r')) and {TV'V" {Q,,r)V"' {Cl' ,r')) . 
In practice however we will consider the Fourier transforms of these objects in 3D, the kurt spectra, which are the power spectra associated with the trispectra, 
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Table 1. Notations 



Power Spectrum (R/F) 
Variance (R/F) 




Eq.lEl) 

Eq-GD 


Bispectrum 


Q& ^^^^ grr'r' 

^1^2^3' ^1^2^3 


Eq.{TD, Eq.(22) 


Skew Spectrum (R/F) 
Skewness (R/F) 


Df^''r"'(ri,r2),Cf^''^"{fci,fc2) 
SV'^"{r^,r2),Sl^'^"{k^M) 


Eq.(33) Eq.JlD 
Eq.(4T), Eq.(42) 


Trispectrum 




Eq.(23) 


Kurt Spectrum (R/F) 
Kurtosis(R/F) 


T-,rr',r"r"' , ^ ^rr',r"r"'/, , n 
(r-l,r-2),C, {ki,k2) 

T-,rr'r",r"' , ^ ^rr'r".r"'/, , n 
(ri,r2),C, {ki,k2) 

„rr',r"r"'. n r,rr'r",r"' , n 

-f^4 {'-l,'-2),-ft:4 (fcl,fc2) 


Eq.l|46l 
Eq.lgO) 
Eq.lHD 


Deconvolution, Mixing Matrix 


Ci,Ci, Mill 1 Cii' 


Eq.{T5), Eq.{T7), Eq.i|39), Eq.(52) 



C;'' " ^ '(A:,fc') andC;^' ' '{k,k').As was the case for the skew spectrum, we not only do harmonic decomposition on the surface of the celestial 



sphere but on as well on the radial direction. For the construction of C 



(rr',r"r"') 



4.2.1 Two-to- Two Kurt Spectrum 



The first of two kurt spectra C\ 



(rr',r"r"') 



{k,k') can be constructed from the harmonic transform of [sTsiT']im.{k) of the quadratic combination of two 



arbitrary spin-weight fields discussed previously in the context of the skew spectrum Eg.iSOt. The resulting kurt spectrum is generic and can be defined for 
any given set of four spin-weight functions defined in 3D. Unlike the skew spectrum which is zero for Gaussian fields, the kurt spectra are non-zero even 
in the absence of any non-Gaussianity, which introduces additional complexity. The Gaussian contribution (also known as the disconnected piece) needs 
to be subtracted out before it is employed for the study of non-Gaussianity. The noise, often assumed Gaussian, can also be subtracted following the same 
technique. It will contribute only to the disconnected part. Later in this section, we will also consider the effect of a mask as we did for the skew spectrum. 

We will use these results to derive expressions forcl^^'^"'^"'>(fc,fc') which leads to: c;^^''^"^"'(fc,fc') = E™[.r.'r']L(fc)^.r",»r"]!m(fc'). 

These power spectra directly probe Tj^''^ {I, k, k'). It compresses all the available information in quadruplets of modes specified by {h, h, Is, U) to a power 

spectrum. The power spectra C^^^ ^ ' (fc, k') and cf^^ ^ ' (fc, k') differ in the way they associate weights to various modes and contain complemen- 
tary information. The reduced trispectrum Ti^i^{ki;L) is defined in terms of {sTj^mi (A:i)s'rj.^„j2 (A;2)s"r"^„^ (A;3)s"'rj"^^ (A:4))c as follows. We have 
added the radial distances associated with each spherical harmonic in the argument with L, which specifies the diagonal formed by the quadruplet of four 

quantum numbers (,ri,„,(fci),,r;^„^(fc2)."r;;„3(fc3).-r;;'„JM>c = Elm(-1)'XV4('=';^^ ^ A Ma U L 
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The final expression depends on the spin indices of various fields as well as on a Kernel F (defined below) which has angular harmonic numbers U and 
radial wave numbers k, k' as its arguments: 
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The kernel F^^''^^ {U, I, k, k') is defined in terms of the reduced trispectrum T^^l^ {ki\ L). 
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We will consider the two components of the trispectrum "snake" and "stars" separately for each of the two kurt specrta. If we follow the algebra, which is 
very similar to what was done previously to derive the skew spectrum we arrive at the following expression for the star component of the three-to-one kurt 
spectrum. The expression reduces to an one dimensional integral as we use the Limber approximation for simplification. 
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The contribution from the star diagram can be expressed in a similar manner. We simply have to replace the /g*^]^^, which is defined in Eq.l|27|l with I^'^l^^ 
defined in Eq.{26j. The additional terms that describe the exchange symmetry of snake terms however will involve the computation of 6j symbols Ea.ll28t 
which poses additional computational complexity. There are two cumulant correlators at four-point level as explained above. 



4.2.2 Three-to-One kurt spectrum 

Following the discussion above we now focus on the other degenerate power spectra associated with the cumulant correlator 
(^rs'r^//r"](f7)s///r"'(f2')>.We win start by defining the all-sky harmonic transform [TT'V"]i^{k) for the cubic field ^rs'r^//r"(n, r) and 
cross-correlate it against the remaining field s"'^im{k')- This is of the same order as {[sV ^iV']{Cl)[^iiV" gii,T"']{Cl')) and contains informa- 
tion about trispectra as well. The compression of the information is done with different weighting for different modes: C, ' {k,k) — 
2lTT5^mR'eal{[s/r^//r"^/"r"'];*„(fc)s/r;„(A:')}- We can now use the definition of the trispectra T^^l^ {L; ki, k2) to express Cl^^'^" {ki, k2) 
in terms of the trispectra. The main difference with the previous spectrum C\ ' {k\,k2) is that it sums over all possible configurations of the 

quadrilateral keeping one of the sides fixed, whereas C\ ' (^1,^2) keeps one of the diagonals fixed but otherwise sums over all possible configuration 

of the quadrilateral. The harmonics of the cubic field can be expressed in terms of the individual harmonics using the following expressions. In the first 
equation we treat the radial direction using real-space expression and next we also do a Fourier transform in the radial direction to obtain a full 3D 
decomposition. 
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Following these expression we can express the three-to-one kurt spectra in terms of the trispectrum. The geometric prefactors that appear in the projected 
(2D) three-to-one power spectra also appear in the 3D expression. The radial harmonics dependence comes through the kernel F^^'^'^: 
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The kernel F'^^'^^ (/^, k, k') is defined in terms of the trispectrum tI^I^ {ki\ L) as follows; 
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The star contribution can be expressed in terms of the kernel I^^l^, using the Limber approximation: 
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For the snake component, we need to replace the kernel J^^^r with -^inakc which also involves a single integration. The evaluation of snake terms are however 
difficult as in addition there will be terms that involve 6j symbols Eg. (1281) because of their permutation symmetry. As was the case with skew spectrum, we 
can collapse these fourth order two-point objects to reduce them to one-point numbers, the kurtosis, which will be a function of both radial harmonics fci , ^2: 



iCr^''^"^"'(fci,fc2) = ^(2/ + l)cf^''^"^"\k^,k2); Kl'"^"'^"'{ki,k2) = ^(2Z + l)cP'^"'^"''(A:i,fo). (51) 
I I 

The corresponding real space (in radial direction) versions K^^ ^ {ri,r2), and K^^ ^ (^1, ''2) can be related using an equation similar to Eq.ll42t. 

We will deal with the mode coupling issues arising from the partial sky coverage next. We will show how to deconvolve the effect of mask. It is 
important to keep in mind that the one point objects such as K^^ ^ (fci, ^2) will have more signal-to-noise and can play important role in studying 
growth of non-Gaussianity along the radial direction. 

A few general comments are in order. At the level of the trispectrum there are two hierarchical amplitudes. If we employ the two sets of kurt spectra 
then the amplitudes can be determined. This is very similar to their use in the CMB where t he kurt spectra can be used to pr ovide independent constraints 
for parameters qnl and tnl that describe the Taylor expansion of the inflationary potential IH1II1999I ; lOkamoto & Hull2002h . Indeed in some hierarchical 
models the snake amplitude Ra and the bispectrum amplitude are related by _Ra = . Similar relation also exists between /at l and rjv l ■ The study of skew 
spectrum can provide direct estimates of the parameter Qg. The estimates from skew spectrum is expected to be more significant statistically due to the higher 
signal-to-noise. An independent estimation using kurt spectra can provide a direct test of various hierarchical ansatz. It is important to realize that the skew 
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spectra as well as the kurt spectra are integrated quantities, i.e. their amplitudes at a specific I depends on the entire range of / values being considered. In 
terms of modelling of these quantities, it means that any successful theoretical prediction will have to correctly model the relevant multi-spectra for the entire 
range of I values being probed. The procedure can be extended to even higher-order. The number of distinct topological diagrams that are needed to build a 
correlation function at a given order will be same as the number of related power spectra, e.g. at fifth order there are three topological amplitudes and three 
multispectra. The procedure outlined above can be extended in such cases as signal to noise of the available data improves. 



4.2.3 The effect of a Mask and Subtraction of Gaussian Contribution 

The partial sky coverage will mean that the measured power spectrum C; (fci , ^2) is not the same as theoretical expectation, but is related as before 

by: c!f^ ^ ' (A:, A:') = Gn'C^f ^ ''^ ' (A;, A:'). In fact it can shown that for arbitrary sky coverage with arbitrary mask the above analysis can be gener- 
alized to arbitrary order of correlation hierarchy. If we consider a correlation function at p + q order, for every possible combination of (p, 5) we will have an 

associated power spectrum. Using the same expression for the mode mixing matrix, we can invert the observed {k\, k^) to C;'''''''(A;i, k^). Hence for an 
arbitrary mask with arbitrary weighting functions the deconvolved set of estimators can be written as: C;'' (A;i, k^) — X^r G'fi^clf'''^ ( '2'i+i ' %+i • 
The matrices G (defined below) depend on the spin indices as well as the power spectrum of the mask. 

visa's", s'" 1 (2' 1) /nl , T \ f ^ ^1 ^ \ f I la ^ \ i |2 ' " "' ^ n 1 I o o 

=i^E^n:^(2^'" + l)(^ (, + ,' + ,") -is + s' + s") )[s"' -s"')\'"^^\' e0,l,±2,3. 



(52) 

However, as pointed out before, if we define kurt spectra c/^'''(A:i, A:2) = {21 + l)cl^'''\ki, ^2) we can use the same mode coupling matrices that 
are used in projection for the purpose of deconvolution. However the 3D treatment introduces a remapping of the radial mode due to the presence of a mask 

C^''(fci,fc2) = E.M„,c;(-'(f^iAi,|^fc2). 

The Gaussian contribution to the trispectrum Q can be written as: 

gllll{ki,k2,ks,ki-L) = (_i)'i+'V(2/i + l){2h + l)C[;-^'(Ai,fc2)CfJ''^"'(A3,fc4)<5Lo5iii,<5!3!4 

H2L + l){-iy'+'^+''S,,iJi^,,Cl^'''\ki,k3)Cl!^'''''\k2,k4) + {2L + l^^^^ (53) 

The various power spectra C;^ (A;i , ks) above include contributions from signal and noise. Next we can compute the Gaussian contributions to CJ 

(rr' 21^"]^"'^ II II 

and C; ' following the same procedure as before just by replacing the trispectrum T^^^^ with its Gaussian counterpart G;J;^ {ki; L). Indeed we will 

have to keep the ordering correct for various U and their counterparts. It is also important to realize that in computing the Gaussian contribution we will 

have to take into account both the signal and the noise Cis (assumed to be Gaussian), i.e, Ci — Cf + Cf . We will denote the Gaussian contributions to the 

(three-to-one) kurt-spectra by t/; ' (A:i, A;2) and (two-to-two) t/j ' (A;i, A;2) respectively: 



lil2l3;L 



g(rr',r"r"')(^^^^^)^ ^ Sa,,,Su,ugllll{l,k^,k2). (54) 

I1I2I3I4, 

For realistic surveys with a mask, the unconnected (Gaussian) contributions to the total kurt spectrum, listed above, can be deconvolved in a manner identical 

to what we have presented before for the connected part of the total trispectrum. The mode mixing matrix for the Gaussian contribution is identical to 

~- (rr'r" r'") ~ (rr' p"p"'j 

what we have introduced in Ea.iSH. From the estimated Vi ' {ri,r2) and ' (''i, ^2) these contributions need to be subtracted out before 

comparing them against the theoretical expectations C|"'" ^ {k, k') = 2l^^f ^ [^IF"' ^ll^j ■ equivalent expression holds for the Gaussian 

contributions that needs to be subtracted: C/;'"'" ''^ ^ {ki_,k2) unAQ^^ ^ {ki,k2). 

We have so far only considered the gravity-induced trispectrum in our discussion. However the kurt spectra for primordial non-Gaussianity can be 
derived in a very similar manner by replacing the gravity-induced trispectrum with a corresponding model for the primordial trispectrum. However it is 
expected that gravity-induced non-Gaussiarrity will dominate the primordial ones at least at lower redshift. 



5 ERROR ANALYSIS 

In the previous section, we have derived the expression for the 3D power spectrum associated with convergence field and indicated how a similar analysis can 
be performed for other spinorial fields. Estimation of these power spectra from noisy data will however will always have to deal with issues such as noise and 
partial sky coverage. An estimator which can deal with such observational constraints was developed for the case of power spectra in Eq.lllSt. It is however 
clear that the estimation of power spectra to be meaningful for any cosmological study we need an approximate handle on the error-bars and their covariance. 
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The error analysis for the PCL estimator that was introduced, can be don e using the formalism used in Munshi et al. 1 2009h whi ch is based on pseudo-Cls 



formalism developed by various authors in the context of CMB data analysis ( lEfstathioul2"004ll2006l.lBrown. Castro & Tavloj2005l) . The contributions to the 
error covariance will have three different components. On large angular scales or small I the error will be dominated mainly by cosmic variance where as the 
high / or smaller angular scale it will mainly be dominated by noise due to the intrinsic ellipticity of galaxies. 

5.1 Power spectrum 

C^'^ {k, k') defines the cross-spectra between two spinorial fields sr(f2) and s'r'(f7). That is Cj ' (fc, k') — ^jqrr X^m RefsP''"* (fc)3'f 'im ( fcPl. It is pos- 
sible t o derive the covariance of estimates under certain simplifying assumptions. The general principles for deriving these results are outlined in lMunshi et al.l 
( l2009h and will not be repeated here, and we quote the results for the ordinary power spectra here. For simplicity we will only consider s = s'. In later sub- 
sections we will also consider covariance matrices for the skew spectrum. The covariance matrix of the estimates is {SC]^'^ {k)SC^,''^ {k')), where SCf'^ (fc) 
are the deviations from the ensemble average {C]^'^ {k)): 

{SC^-^'ik, k')SCl^'{k, k')) ^ E,1?(fc, k') + E,^/^(fc, k') + E,TA(fc, k') 

s -s 



Effik, fc') = ^ { sjc^'^ik, k)Cl'^{k, k)^cY'^'{k', k>)C\r^'{k>, k') + Cf'^' (k, k')C^,'^' {k, fc') I E 

^ ^ ^ lama 

lama ^ ^ 

Sf,r(fc,fc') = ^ E ( 1 -1') {|[^^''^'H!„™. (k', k')[w]t^^J^C^-^{k, k)C^,'^ik, k) + symm. term. 

lama ^ ^ 

Ik^'^^'Hicmc {k, k')[w\i^ma I \lc\-^'{k,k')Cl-^\k,k')y (55) 



+2 

The symm. term, can be constructed by exchanging k and k' as well as F and F'. We have divided the total contribution into three different components. 
The term E,j?(fc, k') is the cosmic variance contribution and depends on the target spectra but is independent of noise. The term E^/^(fc, k') signifies the 
noise contribution and finally E^/^ (fc, fc') is a cross term which gets contributions from both signal and noise. We have assumed that the noise is statistically 
uncorrelated but it varies with pixel-position in the sky; i.e. (F(f2, r)V'{Q,' , r')) = a^'^ {Cl, r, r')52D{^ — H'). 

A detailed modelling of source distribution is required for 3D error estimates. The observational mask w{fl) that we use is completely generic however 
our results uses completeness and orthogonality of spherical harmonics on the cut-sky. This means results will be accurate only for near all-sky coverage. 
The various window functions that we have introduced are constructed from the 3D harmonic transforms such as [a-'^w'^]LM{k) and [aw]LM (k) of maps 
constructed from 3D noise maps and the mask (2D). These window functions are assumed to be much sharper than any variation in the power spectra. 
However such an assumption is unlikely to pose a problem as the weak lensing power spectrum lacks features unlike that of the CMB. The above expression 
is expected to be reasonably accurate at high I regime where the noise dominates. The 3D harmonics that we have used in our derivation based on the following 
definitions: 



[a^''^'{n)w{n)]i„,{k,k') = J [a^'^' {n,k,k')w{Q.)]Yi*m{n)dn 

a^'^ {Q,k,k') = — J kdkji(kr) j k' dk' ji{k'r')a^'^ {d,r,r') — ^ 



2; + 1 



- 2Z + 1 2^ + 1 



(56) 



In the final step we have used the Limber approximation. Similar terms such as [a^'^ (Cl)w{Cl)]im.(k, k') can be dealt with in a similar manner. In practice 
the evaluation of these terms will depend on the redshift distribution of galaxies. 

The deconvolution of the error-covariance matrix can be performed using a similarity transformation. It involves the mode-coupling matrix introduced 
before: 

{5C^'^ {k,k')SCf,'^ {k,k')) = M-^^{5Cl'^ {k,k')SClf {k,k'))M~], A sum over repeated indices is assumed in this equation. The mode-coupling 
matrix introduced here depends on the spins of the relevant fields involved s and s' Eq.dlSb. We will next extend this result to the skew spectrum of arbitrary 
spinorial fields. For higher-order spectra the results are more involved but they can be computed using the same techniques considered here. 

5.2 Skew spectrum 

The expression for the skew spectrum, valid in the high-/ regime using the Limber approximation, is derived in Eq.J35t. The estimator for the skew spectrum 
quoted in Eq.(|38l> depends on cross-correlating an arbitrary product field [sFs'^]im{k) against a third field s"^'i'm{k')- 

In this section we compute the error-covariance under certain simplifying assumptions. We will start from our definition of skew spectrum Eq.l|38t and 
quote the result for the covariance, which can be derived using similar procedure adopted for the derivation of covariance of the ordinary power spectrum. 
The results correspond to three generic fields F, F' and F' with arbitrary spin weight s, s' and s" = s + s' respectively: 

(5Cr''^"(fc, k') k')) « Efi^ik, k') + T.fi^{k, k') + l^^i^ik, k'). (57) 
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The individual terms in terms of noise and signal power spectra are as follows: 

^^2^ (, s + s' -(. + .') )\s" -s" 

x^NN/, , i\ 1 / ' \ I I la I' \ ( < 2 FFxr' i 2 r",r"i* , r 2 rr'.r"/, 7/M2I 

E,,, (fc,fc ) = — ^ Q -(s + s') J s" -s" JV'^'^ \io,mAk)lw a ]i^^^{k) + [wa {k, k )] j 

+ |u''a^^''^^'b„„„(A:)|u>'L„„y^Cf '^"(fc,fc)Cf;"'^"(fc,fc) + |«''a^"''"'L™„ (fc)|i"'|i.m„ \/cf^''^^'(fc, fc')Cf;^''^^'(fc, fc')} (58) 

The error-covariance depends on noise maps for the product field as well as the individual field. The noise in our analysis is not assumed constant and can 
vary with position. The terms such as \w'^a^^'^^ ham^ (k) can also defined using expression similar to Eq.(l56t. It underlines the difficulty associated with 
an accurate error estimation beyond the power spectrum. 

These results are based on various assumptions valid at high -/ regime. However for future, near all-sky surveys for which the harmonic description 
is more appropriate these results can provide a good handle on estimation errors. Alternatively the error-covariance can be computed using Monte-carlo 
simulations. Simulating multiple copies of the observed sky, with all observational details, can be computationally expensive. Hence, often simplifying 
assumptions are employed to compute the covariance. The approach developed here can play a complementary role in cross-checking and validating such 
results. The lower-order covariance such as what we have considered above typically depends on higher-order power spectra. It is customary to quote error 
bars associated with estimated power spectra. However it is important to note that the error-bars for the higher-order spectra such as the skew spectra may 
not be fully informative, as the probability distribution of the skew spectrum for a given I can be skewed. In such cases, the error bars still can give an idea of 
statistical scatter around the estimate. 

The results for deconvolved PCL estimates for the skew spectrum can be computed using a similarity transformation: 

{SCf^'-'^" {k,k')5Cf,^''^" {k,k')) = Ell' GTl (^C'L^''^"('=''=')^C'Lr'^"(fc,fc'))G^,V- ^he mode-coupling matrix G is independent of the radial 
wave vector k but depends on all three spin indices Eg. (139b. Higher-order multi-spectra such as the the skew spectrum are typically more correlated and 
binning may be needed for non-singular inversion. 

It is also possible to compute the cross-covariance of these estimators for power spectrum and the skew spectrum which can be used jointly. These can 
result in tighter cosmological constraints. The results can be derived using the techniques presented above. The coupling matrices are different for different 
spinorial fields. They depend on the spin indices of the constituent spinorial fields. The spinorial fields considered above are however completely generic. If 
we assume that the magentic or B-mode is absent then further simplification can be achieved. 

5.3 Optimal Estimators 

The estimators that we have constructed can be generalised if we optimally weight the harmonics with an inverse variance weight. The generalized two-to-one 
power spectrum Sf^ in this case takes the following form: 

sf^'-^"ik,k') = -±-j2^ Bfiy;,r {k,k,k')BfiT;,y'' ik,k,ky, 
I'l" 

A-^ = ([Cr{k,k)Cr^' (k,k)Cr""' {k',k')]-^ +cyc.peTm.'^ (59) 

This particular result is valid fo r all-sky coverage. The de nominator A is the scatter in the estimator in the Gaussian limit. For partial sky coverage a more 
elaborate treatment in line with lMunshi & HeavensI ( |2009|) is required. In the absence of spherical symmetry due to lack of all-sky coverage or asymmetric 
noise, we will have linear terms in the estimator. The estimator constructed in this way can achieve maximum possible signal-to-noise for a given data 
set. The one-point counterpart for this estimator, denoted as S can be recovered by summing over angular harmonics I, i.e. S^^ {k, k') = E; (2^ + 

1)5^ {k, k'). An interesting point which we note here is that the hierarchical ansatz is factorizable which will allow easy construction of optimal 
weights. In general the expressions for gravity-induced the bi- or trispectra are not factorisable. At bispectrum level we can use the skew spectrum estimator 
to recover the tree amplitude Q3 and similar estimators can be designed for the kurt spectra. The two different kurt spectra will allow independent estimation 
of topological amplitudes Ra and Rt. However it is expected that signal to noise at the level of trispectrum will be low. These optimal estimators can 
be optimized to detect either the gravity-induced non-Gaussianity or different models of primordial non-Gaussianity. They can also be used to forecast 
cross-contamination in a specific estimator from various sources of non-Gaussianity. 

It is worth repeating here that the next generation of weak lensing surveys will have nearly all-sky coverage. This will probe a wide range of angular 
scales. The most commonly used technique for statistical characterization of such surveys is real space analysis using two- or three-point correlation functions. 
One of the advantages of using the higher-order correlation hierarchy is its ability to extract information from a complex survey geometry due to partial sky 
coverage. However, the real space analysis introduces highly correlated measurements for various angular scales. These correlations are more dominant at 
small angular scales where most of the observational information is contained. The alternative is to use the harmonic space representation of the correlation 
functions, e.g. the multi-spectra that are being studied here. Though mathematically equivalent the power spectra or their higher-order generalizations are 
much less used in the context of weak lensing. However the theoretical interpretation of multispectra is much simpler and different harmonics are much less 
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correlated. The main difficulty in harmonic analysis is related to the partial sky coverage. Typically the mask consists of bright stars and saturated spikes 
where no lensing measurements can be performed. The analytical results presented here provide a general analysis of the problem these pose. The results 
relate the convolved and deconvolved power spectra that can be constructed from the higher-order multispectra. The deconvolution process consists of simple 
matrix inversion and can be performed for arbitrary sky coverage. For the case of convergence the matrix representing mode-mode coupling in the presence 
of mask is independent of the order of the multi-spectra being probed. However that is not the case for the case of shear or flexions. The formalism developed 
here also allows for computation of scatter or variance associated with various estimators. 



6 CONCLUSIONS 

It is now well accepted that the next generation of weak lensing surveys will play an important pa rt in further reducing the uncertainty in fundamental 



cosmological parameters, including those that de scribe the evolution of equation of s t ate of dark energy Refregier et alj (120101). They will also be instrumenta l 
in testing various alterna ti ve gravity mode l s (e.g . lHeavens. Kitching & Verdd ( l2007l) : lAmendola. Kunz & Sapond i 20081) : Benvon. Bacon & Kovamal ( l2009h : 



ISchrabback eTai] ( l2009t) : iKilbinger et all ( |2009|) ). The power of weak lensing surveys largely depends on the fact that they can exploit both the angular 
diameter distance and the growth of structure to constrain cosmological parameters. It is therefore very important to develop analytical techniques and 
statistical tools that can fully exploit the potential of future weak lensing surveys. 

Typically without the redshift information the data from weak lensing surveys are analyzed i n projectio n for the entire source distribution. However 
it was found that by binning sources in a few photometric redshift bins the constraints improve l lHiilll999l) . In recent year s a full 3D formalism which 
exploi t s the photomet r ic red shifts of individual sources were developed. Such an approach does not involve any binning; see iHeavensI l [2003l) : ICastro et a3 
l l2005l) :l Heavens et all ( l2006h . Further studies along these line s demonstrate that 3D lensing can provide more po werful and tighter constraints on the dark 
energy equation of state parameter, and on neutrino masses jde Bernardis et al1l2009l : Ijimenez et al. 2O10ll2O10h , as well as testing braneworld and other 
alternative gravity models. These constitute the main science drivers for the future weak lensing surveys. Initial studies in weak lensing focused on two-point 
correlation functions or the power spectrum mainly due to the low signal-to-noise available for higher-order studies fro m most first generation surveys. With 
the av ailability of modern surveys it is useful to include the non-Gaussianity information in the data analysis pipeline jTakada & Jainll2004l : ISemboloni et all 
I2OO9I) that can help to lift some of the degeneracies in estim ation of cosmological parameters involving power spectrum alone. 

In their recent work lMunshi. Coles & Heavenj ( 1201 Ol) have explored the possibility of extending the higher-order statistics of convergence to 3D. The 
main motivation of this work is to gen eralize those results to spinorial ob jects and perform a full 3D analysis for the higher-order statistics. In this sense this is 
also an extension of results derived in iMunshi. Heavens & Coles! ( I2OI0I) which analyzed higher-order statistics of spinorial fields but only in projection (2D). 
The results here are valid for all-sky surveys. It depends on full 3D spherical harmonic decomposition on the surface of the sky as well as along the radial 
directions. Such an approach in analyzing the data from future surveys which will cover a large fraction of the sky. 

The higher-order statistics of convergence k, shear 7± or flexions F and Q depend on accurate modelling of the underlying density contrast 5. Various 
models are used, such as the hierarchical ansatz which we use here, known to be valid in the highly nonlinear regime. However the techniques developed here 
are generic and can also be used in association with other models such as the halo model. 

The higher-order multispectra contain invaluable information. Some of these information is however degenerate because of symmetries associated with 
higher-order correlation functions. It is difficult to estimate the higher-ord er multispectra mode by mode because of the associated scatter involved in such 
estimation especially from a noisy data set. In iMunshi & Heaveiisl ( l2009h . various power spectra (skew spectrum, kurt spectra) were introduced, that are 
associated with a multispectra of a given order and can be estimated in the presence of mask and noise. These spectra carry some of the information contents 
of the multispectra from which they are constructed. In our present study we express the skew spectrum and two degenerate kurt spectra of generic spinorial 
fields in terms of the bi- and trispectrum. This extends earlier results for the convergence (spin-0) field. Extending the previously introduced 3D power 
spectrum C;(fci, ki) to higher-order, we introduce a series of power spectra related to multispectra at each order. We have introduced the 3D skew spectrum 
Cf^ (fci, ^2) associated with the bispectrum of arbitrary triplets of spinorial fields F, F, F". Analogously, at the level of trispectrum we have introduced 
t wo 3D kurt spectra Cf^ ^ (fci , ^2) and Cf^ ^ (fci, k2) for arbitrary choice of spinorial fields. These extends the skew- and kurt spectra defined 
in IMunshi. Coles & Heaveiisl ( I2OI0I) where harmonic decomposition was performed only on the surface of the celestial sphere and a real space analysis was 
performed on the radial direction leading to a mixed representation of skew spectrum Cp'^' (r2, ri) as well as their higher-order counterparts i.e. the two kurt 
spectra C;'^'^^(r2, ri) and CP'"*^' (r2, ri). 

The generic expression for the skew- and kurt spectra involve spherical Bessel functions. We simplified these radial integrals by using the Limber 
approximation, whose accuracy scales as £>( j^). We show that at each order the Limber approximation can reduce the dimensionality of the integrals to unity 
which dramatically reduces computational cost. Both the Limber approximation and the hierarchical ansatz are accurate at smaller scales and their joint use 
can help us to compute the skew- and kurt- spectra very efficiently with reasonable accuracy, but the method can accommodate different models for nonlinear 
clustering. 

We also present analytical results for dealing with a mask, via a pseudo-C; approach, encapsulated in a mode-mixing matrix. The estimation of unbiased 
skew- or kurt spectra are done by simple inversion of the mixing matrix M, which depends on the spins associated with the spinorial fields. Some regularisation 
will normally be required. The presence of an observational mask typically only induces mode-mixing on the celestial sphere and not on the radial direction. 
We have also showed how our formalism presented here can be used also for the computation of scatter under certain simplifying assumptions in the presence 
of an observational mask, and we have identified individual terms that correspond to contributions from noise, partial sky coverage (cosmic variance) and 
cross terms. 

The results presented here will be relevant for the study of cosmic magnification studies in 3D as well as in many other contexts where integrated radial 
information is used. The estimators for skew or kurt spectra that we have described here can be improved by inverse variance weighting of 3D harmonics. 
Finally, to summarize: 

• We have studied higher-order multispectra in the context of 3D weak lensing surveys. 

• We use a full 3D Fourier decomposition which employ spin-weight spherical harmonics. 
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• Our generic results are valid for arbitrary 3D spinorial objects. 

• The results are relevant for convergence k, magnification ji, shear 7± as well as flexions T and Q or an arbitrary scalar tracer field $. 

• In our analysis we define power spectra Ciiki, k2) that are related to the bispectrum (skew spectra) and to the trispectrum (kurt spectra). 

• We provide both all-sky exact results and corresponding approximate results using the Limber approximation. 

• Use of Limber's approximation reduces multidimensional integrations along the radial direction to one-dimensional integrals. 

• We show how the multi-spectra can be recovered from a masked sky in the presence of noise, and show how the presence of masks mixes modes not 
only on the surface of the sky but in the radial direction. 

• The modelling was done using the hierarchical ansatz but the formalism can work with any input underlying density multispectra. 

• Under certain simplifying approximations, we also obtain expressions for the covariance of our power spectra and skew spectra estimators. 

• We outline how inverse variance weights can be introduced and optimal estimators can be defined for the detection of a specific type of non-Gaussianity. 

• The formalism can be relevant in many other contexts where line-of-sight integrations of non-Gaussianities are performed or in studies involving cross- 
spectra or mixed-bispectra. 

In this paper we have ignored many observational complexities for simplicity, such as that in a realistic survey the lensing potential can only be sampled 
at the discrete positions of galaxies, and the average number of source galaxies will decline with redshift. We also ignore photometric redshift errors. 
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APPENDIX A: USEFUL MATHEMATICAL RELATIONS 
Al Spherical Bessel Functions 

The orthogonality relationship for the spherical Bessel functions is given by the following expression: 
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k^3i{kri)ji{kr2)dk - 



2rl 



SiD{ri — r2). 



The extended Limber approximation is also implemented through the following approximate relation lLo Verde & AfshordillioO^ : 
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Thus, for high /, the spherical Bessel functions can be replaced by a Dirac delta function Sid'- 



A2 3j Symbols 

The following properties of 3j symbols were used to simplify various expressions. 
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